function out=varianceTheoretical(funcmod,parVec,solveOptions,additionalOptions)
% Compute theoretical variance and std of observables and states 
%
% function out=varianceTheoretical(funcmod,parVec,solveOptions,additionalOptions)
%
% Output
% out.varStates     =lyapunov_symm(GG, RR*( SDX'*SDX)*( RR') ); 
% out.varObservables=ZZ*varStates*ZZ'; 
% out.stdStates     =sqrt(diag(varStates)); 
% out.stdObservables=sqrt(diag(varObservables)); 
% out.ACZeroStates  = out.varStates./( out.stdStates*out.stdStates'  ); 
% out.ACZeroObservables= out.varObservables./( out.stdObservables*out.stdObservables' ); 
% AJ Dec 2011 
% 
% 1. Solve Model 
[GG,RR,~,eu,SDX,ZZ]=...
    feval(funcmod,parVec,solveOptions,additionalOptions); 
if ~isequal(eu,[1;1]) 
    error('Model does not have a solution') 
end 

% 2. Compute Theoretical Variance 
out.varStates     =lyapunov_symm(GG, RR*( SDX'*SDX)*( RR') ); 
out.varObservables=ZZ*out.varStates*ZZ'; 
out.stdStates     =sqrt(diag(out.varStates)); 
out.stdObservables=sqrt(diag(out.varObservables)); 
out.ACZeroStates  = out.varStates./( out.stdStates*out.stdStates'  ); 
out.ACZeroObservables= out.varObservables./( out.stdObservables*out.stdObservables' ); 